Simulated Annealing with Tsallis Weights 
A Numerical Comparison 



Ulrich H.E. HansmannF] 

Department of Theoretical Studies 
Institute for Molecular Science 
Okazaki, Aichi 444 > Japan 



ABSTRACT 

We discuss the use of Tsallis generalized mechanics in simulated annealing algorithms. 
For a small peptide it is shown that older implementations are not more effective than 
regular simulated annealing in finding ground state configurations. We propose a new 
implementation which leads to an improvement over regular simulated annealing. 
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Simulated annealing |lj has become an often used tool to tackle hard optimization 
problems in various fields of science. Its underlying idea of modeling the crystal grow 
process in nature is easy to understand and simple to implement. Any Monte Carlo or 
molecular dynamics technique can be converted into a simulated annealing algorithm by 
allowing the temperature to decrease gradually. However, the method is not without 
problems. The performance of simulated annealing depends crucially on the annealing 
schedule. It could be shown that convergence to the global minimum can be secured for a 
logarithmic annealing schedule |2j , but this is of little use in applications of the method. 
Constraints in available computer time enforce the choice of faster annealing schedules 
where success is no longer guaranteed. Similar to the growing of real crystals which hardly 
ever can be archived by a simple cooling process, elaborated and often system specific 
annealing schedules are frequently necessary to obtain the global minimum in the CPU 
time available. Hence, ever since the seminal article by Kirkpatrick et al. [|IJ, attempts 
were made, to improve the performance of simulated annealing in practical applications, 
see for instance Ref. ||. Two of the more recent attempts @, [|| are inspired by Tsallis 
generalized mechanics || [7j . We compare the performance of both implementations with 
regular simulated annealing by taking an energy function for the protein folding problem as 
an example. We describe then a new application of Tsallis weights to simulated annealing 
which leads for our model to an improvement over regular simulated annealing. 

In the Tsallis formalism ||, a generalized mechanics is constructed by maximizing a 
generalized entropy 



S = -kY,Pi^Pi 




with the constraints 
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Here, q is a real number. A generalized probability distribution 
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follows and the average of an observable O can be defined by 



<o>=5>?0 4 . 
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It is evident that for q — > 1 the generalized distribution of Eq. [3] tends toward the Gibbs- 
Boltzmann distribution and therefore regular statistical mechanics is recovered in this 
limit. 

The important feature of Tsallis generalized statistic for optimization problems is that 
the probability of states does no longer decrease exponentially with energy but according 
to a power law where the exponent is determined by the free parameter q (see Eq. [|). This 
observation inspired a generalized simulated annealing algorithm ||] where the acceptance 
probability 



was introduced. Here, AE is the change in energy. Again, for q — > 1, the acceptance 
probability of canonical simulated annealing is recovered. The algorithm was employed 
to find close to optimal solutions to the traveling salesman problem and it was claimed 
that it is faster than classical simulated annealing and has optimal performance for large 
negative parameters q [f|]. However, the above algorithm does not obey detailed balance 
and therefore convergence to an equilibrium distribution is in general not guaranteed. For 
this reason, it was recently proposed to utilizes the acceptance probability || 



with lirciT^o q{T) = 1. Detailed balance is obtained by the algorithm and convergence to 
the generalized distribution of Eq. [3] is guaranteed for each temperature T and parameter 
q. Since q — > 1 as T — > 0, this generalized simulated annealing algorithm tends in the 
same way as regular simulated annealing to a steepest descent at low temperature. For 
Tsallis parameters q{T) > 1 the distribution of Eq. |3| has a tail to higher energies and 
the probability to cross energy barriers and to escape local minima is therefore increased 
by the above weights. The new algorithm was applied in Ref. || to the conformational 
optimization of the 48 atom tetraalanine peptide using molecular dynamics and hybrid 
MD-MC methods in the CHARMM force field [§. Both temperature and the Tsallis 
parameter q were exponentially decreased, with a start value of q = 2 for the Tsallis 
parameter. Results better than conventional simulated annealing were reported. 



i -i 



(5) 



p(AE) = min 1, [1 - (1 - q)/3AE] 
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However, a disadvantage of the above algorithm is that it requires the careful tuning of 
additional free parameters. Not only a suitable annealing schedule in temperature T has 
to be chosen, but also one in the Tsallis parameter q(T). Furthermore, the performance 
of this algorithm depends also on the choice of the zero in potential energy (which was 
ignored by the authors of Ref. ||). This is because the transition probabilities of Eq. |^ 
are not invariant under a shift in energy. For instance, even for q > 1 (the case q < 1 
can lead to complex probabilities) the weights of Eq. ^ can become negative for negative 
values of the energies. On the other hand, acceptance of a configuration will no longer 
depend on the energy of the configuration if all energies are shifted by a large enough 
positive number. Hence, while classical simulated annealing has already the problem of 
finding an optimal annealing schedule for the temperature, the above algorithm requires 
in addition a careful tuning of the Tsallis parameter q and of the energy scale, making its 
applicability highly model dependent. 

Here we show how the above problems can be alleviated. Use of Tsallis weights in 
the course of a simulated annealing simulation is motivated by the fact that the resulting 
probability distribution has a tail to higher energies for q > 1, enhancing in this way the 
probability to cross barriers and escape local minima. While negligible at high tempera- 
tures this feature becomes important at low temperatures where a canonical weight would 
make it difficult to escape local minimas. We are therefore mainly interested in the use of 
of Tsallis weights for low temperatures. It is obvious that the Tsallis distribution at low 
temperatures should not be dominated by the tail to higher energies, but still be centered 
around the energy where the canonical distribution has its maximum. For otherwise, low 
energy (temperature) states will not be sampled sufficiently To find the parameter q > 1 
which yields to an optimal distribution let us first write the Tsallis weights as 

q 

w (E) = [l-(l-q)l3(E-E )]l-<l , (7) 

where Eq is the (in general unknown) ground state energy. This is equivalent to chos- 
ing an energy scale where all energies are positive with the zero for the ground state. 
In this way we ensure that the weights are always positive. The Tsallis weights will 
be a good approximation of the Boltzmann weights Wb(E) = exp(—j3(E — E )) for 
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(1 — q)(3(E — Eq) « 1 . To ensure that simulations are able to escape from energy local 
minima, the weights should start deviating from the exponentially damped Boltzmann 
weights at energies near its mean value. This is because at low temperatures there are 
only small fluctuations of energy around its mean. We may thus set in Eq. ^ 

-(l-q)P(<E> - E ) = ± (8) 

The mean value of energy is given at low temperatures by the harmonic approximation: 

< E > -E « ^kT = %, (9) 

where np is the degree of freedom of our molecule. Hence, for low temperatures, Eq. |^ 
can be written, as 

-<i- 9 >=H, do) 

which leads to an optimal Tsallis parameter 

<? = ! + —• (11) 
n F 

Hence, we propose a generalized simulated annealing algorithm where configurations 
are weighted with 

_Q 

w( 



,(E) = [l-(l- q ){3(E-E )]±-<l . (12) 

The Tsallis parameter q is set to q = 1 + 1/np. Our weights require knowledge of the 
ground state energy. However, in general E is not known. We therefore approximate 
E in the course of a simulated annealing simulation by Eq = E min — c where E min is 
the lowest energy ever encountered in the simulation and c a small number. Eq is reset 
every time a new value for E m i n is found. Changing the value of Eq is a disturbance of 
the Markov chain and while we expect the disturbance to be small, we clearly cannot use 
our algorithm to calculate thermodynamic averages. However, due to finite stepsize of 
the temperature annealing we can anyway not assume convergence against an equilibrium 
distribution. As with regular simulated annealing, our method is valid only as a global 
optimization method. 

We have tested the various simulated annealing algorithms for the protein folding 
problem, a long-standing problem in biophysics with rough energy landscape. Here, Met- 
enkephalin has become a often used model to examine new algorithms. Met-enkephalin 



has the amino acid sequence Tyr-Gly-Gly-Phe-Met. The potential energy function E to t 
that we used is given by the sum of electrostatic term Ec, Lennard- Jones term Elj, and 
hydrogen-bond term E^b for all pairs of atoms in the peptide together with the torsion 
term E tors for all torsion angles. The parameters for the energy function were adopted 
from ECEPP/2.[y] Fixing the peptide bond angles u to 180° leaves us with 19 torsion 



angles as degree of freedom. The computer code KONF90 [10| was used. 



As in earlier work on Met-enkephalin [] 1 we performed for each of the different algo- 



rithms 20 runs of 50,000 sweeps. Each run started from completely random configuration 
and each angle is updated once in a sweep. The temperature was lowered exponentially 
according to 

T = TstY" 1 (13) 

for the ith sweep and 

7 = (T FI /T ST )^ , (14) 

where T$t is the start temperature and Tpj is the final temperature. One of the quanti- 
ties we monitored to evaluate the performance of the various algorithms was the average 
< E Low > (taken over all 20 runs) of the lowest energies E Low obtained in each single run. 
The other quantity was the number ng of ground-state configurations found in the 20 



independent runs. In Ref. [|L2| it was shown that with the energy function KONF90, con- 
formations of energy less than —11.0 kcal/mol have essentially the same three-dimensional 
structure. Hence, we consider any conformation with E < —11.0 kcal/mol as the ground- 
state configuration. 

Let us present now our results. In Tab. 1 we show typical results for the first general- 
ized simulated annealing algorithm which was proposed in Ref. || and uses the acceptance 
probability of Eq. |5[ We chose as start temperature T$t = 1000 K and tried two values 
for the final temperature: Tp = 50 K and Tp — 1 K. These temperatures were also 



chosen by us in Ref. [JT/TJ from which we have also taken the results for regular simulated 
annealing (indicated by q = 1 in Tab. 1). We did not observe an improvement over 
regular simulated annealing for any choice of the Tsallis parameter q in this generalized 
simulated annealing algorithm. In the range < q < 1.25 the performance is comparable 
with canonical simulated annealing. Both uq and < El ow > vary little in this range. The 
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performance of the algorithm detoriates quickly outside of this range. This is especially 
true for large negative values of q which were presented as the optimal choice for q in 
Ref. ||. The probability of finding ground state configurations becomes small and the 
average < E Low > of lowest energies found in each single run not only increases, but also 
the standard deviation of < Ep ow > indicating that the obtained low energies strongly 
depend on the initial configuration and the algorithm reduce in this case to a mere quench- 
ing. Hence, for our system, application of this generalized simulated annealing algorithm 
seems to bring no improvement over regular simulated annealing. 

In Tab. 2 we show our results for the second method, which uses the acceptance 
probability of Eq. |^ and was proposed in Ref. ||. Following the authors of Ref. || 
we chose the start value q = 2 and decreased both temperature and Tsallis parameter 
exponentially such that for the final temperature Tpj we have q(Tpj) = 1. Again we 
tried two values for the final temperature: Tp = 50 K and Tp = 1 K. Since Tsallis 
distributions have a tail to high energies for q > 1, we expected that it would be possible 
to choose lower start temperatures than for regular simulated annealing. Hence, we tried 
not only T$t = 1000 K, but also T$t — 500 K and 300 K. In each case we found a poorer 
performance than for regular simulated annealing. This is in contradiction to the results 
of Ref. [Q, where the authors found a significant improvement over canonical simulated 
annealing. We remark that for the calculation of weights (see Eq. ||) in the simulations, 
we shifted the energies by the ground-state energy for Met-enkephalin (as known from 
previous work [[□]]) to ensure positive weights, otherwise even poorer results were obtained. 
Hence, the poor performance of this algorithm in simulations of our peptide is not only 
due to a poor choice of the energy scale, but also to an imperfect annealing schedule of 
the Tsallis parameter q(T). We conclude that the performance of this algorithm is highly 
model dependent and requires carefully tuning in the annealing of q and the choice of the 
energy scale. This is a severe limitation for practical applications. 

Finally, Tab. 3 shows the results for our implementation of Tsallis weights in simulated 
annealing algorithms using the acceptance probability of Eq. |12[ The main difference to 
the previous algorithm is that the Tsallis parameter q is not free, but set to an optimal 
value of qp = 1 + l/n F = 1 + 1/19 ~ 1.053. In addition, our procedure guarantees that 
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the weights are always positive by self-tuning the estimate of the ground state energy Eq 
in the course of an annealing run. Eq is reset every time to Eq = E m i n — 1 kcal / mol when 
a new configuration with lower energy E min than any previous configuration is found. We 
found for both canonical and generalized simulated annealing an optimal performance for 
start temperature T$t — 500 K and final temperature Tp = 50 K. With this temperature 
annealing schedule the ground state configuration was found 8 out of 20 runs for regular 
simulated annealing and 12 out of 20 runs for generalized simulated annealing. This 
is a modest improvement of the new algorithm over the canonical simulated annealing. 
The improvement can also be seen in the estimate for < E^ ow > which is 0.6 kcal/mol 
lower for the new algorithm and has a smaller standard deviation than regular simulated 
annealing. We further notice that the new generalized ensemble algorithm allows to start 
the temperature annealing at lower temperatures. While regular simulated annealing 
works best with start temperatures over 500 K, the performance of the new algorithm 
depends only little on the start temperature and rather favors T$t < 500 K. This 
follows from the form of the Tsallis distributions which have a tail to high energies for 
q > 1. Equilibrization at lower temperatures is therefore enhanced. We remark that the 
improvement of the new method over regular simulated annealing still does not lead to the 
performance reported for the generalized ensembles algorithms in Refs. [II], O] . However, 
the new algorithm is much easier to implement. Since the Tsallis parameter q is constant 
in our algorithm, only minor modifications are required in existing simulated annealing 
programs to accommodate the new technique. Unlike in the algorithm of Ref. || the 
improvement over regular simulated annealing is gained without the need of determining 
optimal annealing schedules for additional parameters. 

Let us summarize our results. We have performed Monte Carlo simulations of Met- 
enkephalin using Tsallis generalized Mechanics. We discussed older proposals for the 
use of Tsallis weights in simulated annealing algorithms and showed how to overcome 
their shortcomings. Our algorithm is easy to implement in existing simulated annealing 
programs and does not require tuning of annealing schedules in additional variables. For 
the case of a simple peptide we have demonstrated that our new technique offers an easy 
way to improve the performance of simulated annealing. 
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Table 1: 

Results for simulated annealing simulations using Tsallis weights as defined in Ref. J3J. 
For each value of q 20 independent runs of 50,000 sweeps were performed. The start 
temperature was T$t = 1000 K. q — 1 indicates results from regular simulated annealing 



which were taken from previous work in Ref. JTT]. n G is the number of runs where a ground 
state configuration was found. < E^ ow > is the average over 20 runs of the lowest energies 
obtained in each run. The standard deviation of this quantity is given in parentheses. 





T F 


= 50 K 


T F = 1 K 


1 


n G 


< E Low > 


n G 


< E Low > 


2.00 


0/20 


-6.0 (0.9) 


2/20 


-9.2 (1.3) 


1.75 


0/20 


-7.5 (1.0) 


1/20 


-9.5 (1.2) 


1.50 


3/20 


-9.6 (1.1) 


1/20 


-8.7 (1.3) 


1.25 


5/20 


-9.5 (1.2) 


4/20 


-9.7 (1.4) 


1.00 


6/20 


-10.0 (1.3) 


8/20 


-10.0 (2.2) 


0.75 


6/20 


-10.2 (1.3) 


7/20 


-9.9 (1.9) 


0.50 


5/20 


-10.2 (1.3) 


6/20 


-9.9 (1.6) 


0.25 


5/20 


-10.2 (1.3) 


8/20 


-10.0 (1.7) 


0.00 


8/20 


-10.4 (1.3) 


2/20 


-9.0 (1.8) 


-0.50 


5/20 


-9.3 (1.8) 


5/20 


-9.1 (1.9) 


-1.00 


4/20 


-9.5 (1.8) 


5/20 


-9.6 (1.9) 


-2.00 


1/20 


-8.6 (1.8) 


1/20 


-8.1 (1.9) 
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Table 2: 

Results for simulated annealing simulations using Tsallis weights as defined in Ref. ||. 
For each annealing schedule characterized by the choice of start temperature T$t and 
final temperature Tp 20 independent runs of 50,000 sweeps were performed. The results 
are compared with that of regular simulated annealing runs, n G is the number of runs 
where a ground state configuration was found. < E Low > is the average over 20 runs of 
the lowest energies obtained in each run. The standard deviation of this quantity is given 
in parentheses. 



Tst/K 


T FI /K 


Reg 


ular Simulated Annealing 


Simulated Annealing Version of Ref. 5 






n G 


< El ow > 


n G 


< El ow > 


1000 


50 


6 


-10.0 (1.3) 





-7.9 (1.8) 


1000 


1 


8 


-10.0 (2.2) 





-7.9 (1.3) 


500 


50 


8 


-10.5 (1.3) 


3 


-8.3 (1.8) 


500 


1 


2 


-9.3 (1.3) 


1 


-8.1 (1.6) 


300 


50 


1 


-9.8(1.2) 


1 


-8.3 (1.4) 


300 


1 


3 


-9.6(1.4) 


5 


-8.9 (2.0) 
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Table 3: 



Results for simulated annealing simulations using Tsallis weights as defined in Eq. 12. For 
each annealing schedule characterized by the choice of start temperature T$t and final 
temperature T F , 20 independent runs of 50,000 sweeps were performed. The results are 
compared with that of regular simulated annealing runs, n G is the number of runs where 
a ground state configuration was found. < E Low > is the average over 20 runs of the 
lowest energies obtained in each run. The standard deviation of this quantity is given in 
parentheses. 



Tst/K 


T FI /K 


Reg 


ular Simulated Annealing 


New Simulated Annealing Version 






n G 


< E Low > 


n G 


< E Low > 


1000 


50 


6 


-10.0 (1.3) 


7 


-10.7 (0.9) 


1000 


1 


8 


-10.0 (2.2) 


7 


-10.7 (1.3) 


500 


50 


8 


-10.5 (1.3) 


12 


-11.1 (0.9)) 


500 


1 


2 


-9.3 (1.3) 


11 


-10.9 (1.3)) 


300 


50 


5 


-10.1 (1.3) 


13 


-11.0 (0.9) 


300 


1 


3 


-9.6 (1.4) 


11 


-11.0 (1.1) 
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